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Discrete scale invariance (DSI) has recently been documented in time-to-failure rupture, earth- 
quake processes and financial crashes, in the fractal geometry of growth processes and in random 
systems. The main signature of DSI is the presence of log-periodic oscillations correcting the usual 
power laws, corresponding to complex exponents. Log-periodic structures are important because 
they reveal the presence of preferred scaling ratios of the underlying physical processes. Here, we 
present new evidence of log-periodicity overlaying the leading power law behavior of probability 
density distributions of affine random maps with parametric noise. The log-periodicity is due to 
intermittent amplifying multiplicative events. We quantify precisely the progressive smoothing of 
the log-periodic structures as the randomness increases and find a large robustness. Our results 
provide useful markers for the search of log-periodicity in numerical and experimental data. 
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I. INTRODUCTION 

Complex critical exponents and complex fractal dimen- 
sions have until recently been discussed only for hierar- 
chical systems, be they man-made or naturally oc- 
curring as in the mammalian bronchial tree ]9|JlO| . These 
hierarchical systems are characterized by discrete scale 
invariance (DSI), a notion qualitatively similar to the 
concept of "lacunarity" . A signature of this DSI is the 
presence of log-periodic oscillations correcting the usual 
power laws, corresponding to complex exponents. 

Recently, their occurrence in irreversible rupture 
|l6| and growth processes 0,|l||] as well as prior to 
financial crashes |l9|,p0[ have been documented. It has 
been suggested []2 lfj that complex exponents are rather 
common and should be looked for generically in any 
model whose critical properties are described by an un- 
derlying non-unitary field theory. This excludes the 
usual homogeneous spin systems in which the renormal- 
ization flow is a gradient |^2| . This includes models 
with non-local properties such as percolation and an- 
imals 
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polymers and their generalizations, models 
of irreversible growth processes such as rupture [ pd| fl6| , 
Diffusion-Limited-Aggregation (DLA) |l7|]is}|, a nd mod- 
els with quenched disorder like spin-glasses p3|^7|. See 
p8| for a review. 
Three outstanding problems remain : 

• do we know all the physical mechanisms that can 
produce complex critical exponents? 



• how strong are the log-periodic structures and how 
robust are they with respect to noise and disorder? 

• does there exist a smooth invariant probability dis- 
tribution (having a density), or is it discrete? 

With respect to the first question, six situations have 
been discussed : 

1. the presence of a built-in geometrical hierarchy 

§10]; 

2. the diffusion in anisotropic quenched random lat- 
tices in which the hierarchy is constructed dynami- 
cally due to the probabilistic encounters with traps 

{§; 

3. intermittent amplification processes [ ^0[ ; 

4. cascades of ultra-violet instabilities as in rupture 
and growth processes p7|Jl8| ; 

5. non-local geometry J2lj ; 

6. quenched disordered systems [p3| 

In regards to the second question, log-periodic oscil- 
lations of spin systems on a fractal amounts to exceed- 
ingly small effects, typically of the order of 10~ 5 in rel- 
ative value (^||. In contrast, it is still not fully un- 
derstood why log-periodic structures seem to be many 
times stronger, of the order of 10% or so, in rupture and 
growth processes. In addition, log-periodicity implies a 
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prefered scaling ratio which, in nature, should be largely 
perturbed by disorder. Theoretical estimates of the effect 
of disorder on the log-periodic corrections indicate that 
they should be generally robust |Hj . An important prac- 
tical question is how much disorder or noise will make the 
log-periodic corrections too small to be observed. Ensem- 
ble averaging is also an issue as finite-size effects cause 
significant variations in the phase of the log-periodic os- 
cillations. Averaging may cause them to disappear. This 
was observed in DLA clusters [|l7) in which single clus- 
ter analysis uncovered the log-periodic structures while 
averaging procedures destroyed them. 

In order to address these questions on the effect of 
disorder, we study a simple, positive, random map with 
parametric noise , 



X t+1 = a t X t + b t , with a t ,b t > 



(1) 



The growth rate at and the additional term b t are as- 
sumed to be pairs of positive identically distributed ran- 
dom values with the joint distribution function V a ,b- In 
most of the cases treated below, we assume that at and 
bt are independent, which yields V a b — PaPb- 
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FIG. 1. 



Xt history for at with a two point distribution 
(section IV) at a = 2, £ = 2, p — 0.95, and b t = 1. 



It may seem that the linear model ([!]) is so simple that 
it does not require a careful mathematical investigation. 
This is however not the the rather extensive 

mathematical analysis of the problem in |33| indicates. 
We will show a very unusual behavior of solutions of the 
difference equation ([!]). It is known |33) that, provided 
some regularity assumptions hold and when the average 
rate of growth (In a*) is negative, then the time series X t 
is stationary. Furthermore X t is characterized statisti- 
cally by a probability distribution function with a power 
law tail : 



Vx{x) ^ x- {1+ ^ 
when the equation for fj,, 



(2) 



(3) 



has a positive solution |31j |3 

The power law distribution function stems from an 
intermittent and transient amplification occurring when 
several successive at are larger than 1. Its origin is thus 
in the class of intermittent amplifications ]3Cfl and in- 
termittent trapping 29 mechanisms. We may therefore 
expect complex valued \i exponents. This should lead to 
the occurrence of detectable log-periodic corrections in 
the leading simple power law behavior. 

We are here concerned with the continuity or discrete- 
ness of this distribution function and with the strength 
and detection of the potential log-periodic corrections, es- 
pecially as a function of the distributions V a and Vb- In- 
tuitively, the broader these distributions are, the weaker 
we expect the log-periodic corrections to be since a log- 
periodicity is the signature of a favored scaling ratio. 
This preference must disappear as the disorder increases. 
We aim to carefully quantify this scenario, for the benefit 
of future analysis of log-periodicity. 
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FIG. 2. Cumulative distribution of the 10 4 largest iter- 
ates among 10 8 realized for at with a two point distribution 
(section IV) at a — 2, p = 0.95, and b t = 1. 



In the next two sections, we recall useful information, 
discuss a connection with products of random matrices 
and iterated function systems and review the so-called 
transition operator approach determining the probabil- 
ity density function (pdf) P(X). We then discuss the 
case in which at take only two values 1/a and with 
probability p and 1 — p respectively, where a > 1 and 
£ > 0, first in the case of a fixed 6=1 and then with 
increasingly widening Vb distributions. We then analyze 
the case in which a t is broadly distributed and discuss 
the detection criteria for the log-periodic corrections. 

In addition to the present focus as a paradigm of sys- 
tem exhibiting complex exponents, this random map (|l|) 
has been introduced in various contexts, for instance 
in the physical modelling of ID disordered systems |3l]] 
and the statistical representation of financial time series 
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p2| ■ The variable X t is known in probability theory as a 
Kesten variable Q . The map (|l]) describes for instance 
the time evolution of a fish population X t with a* de- 
pending on the rate of reproduction and on the depletion 
rate due to fishing as well as environmental conditions 
and bt describing the input due to restocking from an 
external source such as a fish hatchery, or from migra- 
tion from adjoining reservoirs |3q| . The random map 
(|l|) can also be applied to other problems of population 
dynamics, epidemics, investment portfolio growth, and 
immigration across national borders |3q]. Variations of 
this model have recently been proposed for the analysis of 
crop control in the presence of weed infestion |$7| . Mod- 
els of enonomic evolution typically involve a system of 
affine coupled equations of the type (fO) below, which are 
multi-dimensional generalization of For instance, the 
economic model of Keynes in its simplest form links con- 
sumption, investment and production in a linear affine 
system of deterministic equations. The system ( |Io| ) cor- 
responds to a generalization in which the coefficients of 
the auto-regression are allowed to fluctuate in time to 
account for uncertainty. More generally, models used in 
econometrics |Q are very similar to (|f|) and (JfTJ) below, 
even if they usually assume constant coefficients. 

It is probably true that ([!]) is one of the simplest linear 
stochastic equation that can provide an alternative mod- 
eling strategy for describing complex time series. We 
note that a nonlinear version with a quadratic nonlincar- 
ity (corresponding to the logistic equation with random 
multiplicative noise) has recently been shown to lead to 
a new type of crisis, in that there is a sudden qualita- 
tive change in the chaotic dynamical behavior induced 
by variations of the parameters J3§]. We do not discuss 
these properties but restrict our considerations to the 
affine random map (El). 



II. RESULTS ON THE KESTEN AFFINE 
RANDOM MAP 

A. Formal solution 

The formal solution of (Q) for N > 1 reads 



N-l 



N-l N-l 



Xt+N = ( n a t+l )x t + j2 b t+i n • 

1=0 m=l+l 



>t-\-m ; 



(4) 



2=0 



where we define Jlm=]v a t+m = 1 f° r t ne special value 
I = N — 1. It is clear that the flto 1 at + l multipliers of 
(|J) control the X t dynamics. Thus, X t +N diverges (resp. 
remains bounded) if the average logarithmic growth fac- 
tor (In at) is positive (resp. negative). Here, we focus our 
attention on the case 



In this regime, we notice the role of bt which provides a 
reinjection mechanism allowing X t to fluctuate with- 
out converging to zero, as it would if bt vanished. 



B. Product of random matrices 

The map ([!]) can be written as a product of random 
2x2 matrices : 



Xt+i 
1 





( a t 






)- 


v o 







(6) 



By Furstenberg's theorem, the norm 1 1 V t \ | (~ X t for large 
X t ) of the t-th vector 



V t = 



X, 
1 



\\Vt\\ = \\V \\e^ 



(7) 



(8) 



where Ai is the largest Lyapunov exponent of the product 
of the random matrices. The 2x2 matrices are triangular 
and thus 



Ai = max{(lnat), 0} . 



(9) 



(In a*) < 



(5) 



We recover the exponential growth regime of X t for 
(In a t ) > 0. In the reverse case, (lna t ) < 0, the Lyapunov 
exponent is zero, which corresponds to the marginal case 
between exponential growth and exponential decay. This 
is the regime where one usually encounters power law 
behavior, for instance in power law sensitivity to initial 
conditions in dynamical systems at the onset of chaos 
©■ _ 

It is worth noticing that this zero Lyapunov expo- 
nent is different from the directly measured Lyapunov 
exponent of ([!]). Indeed, the solution (|4|) shows that 
a perturbation SX t at time t gives an error SX t +N = 
SX t Ufji 1 <k+l ~ e( JV - 1 )< lna >. This corresponds to a 
negative Lyapunov exponent for the case studied here 
(H), equal to (In a). This would lead one to conclude 
that the dynamics is trivial. Generally speaking, the 
widespread opinion in the physical community that no- 
tions like chaos, etc., have some strict correspondence to 
the positivity of Lyapunov exponents is not quite cor- 
rect. See, for example, the detailed discusssion of non- 
chaotic dynamical systems with positive Lyapunov ex- 
ponents and vice versa, and futher references in |^3|,Q. 
Here, the usual calculation of the Lyapunov exponent is 
not sensitive to the "reinjection" mechanism brought by 
the bt term. By construction, the matrix formulation (|^) 
takes this effect into account. The resulting vanishing 
Lyapunov exponent alerts us to the possibility of com- 
plex behavior. 
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C. Iterated Function System 

We also mention the relationship with Iterated Func- 
tion Systems (IFS), which are defined as follows [5q| . One 
first defines an affine transformation W from R u to R D : 



W[x] = Ax + b 



(10) 



where A is a D x D matrix and b a vector in R D . An 
affine transformation is contractive if there exists a Lip- 
schitz constant s < 1 such that 



\W{x]-W[y}\ <s|x-y| 



(11) 



An IFS consists of N affine transformations Wi and a set 
of probabilities Pi > with J2i=i Pi = 1- Starting with a 
given set of points, the IFS code consists in applying to it 
an infinite sequence of transformations, each of them be- 
ing chosen with its corresponding probability. In general, 
IFS codes satisfy the average contractive condition : 



< 1 



(12) 



Taking D — 1, we see that (|l0|) is the same as (|Tj), where 
N is the number of different values taken by at (suppose 
for simplicity that bt is constant) with their respective 
probabilities pi. In other words, the affine random map 
(|l|) is a one-dimensional IFS. Then, the Lipschitz con- 
stant Si is equal to the i-th value a, that at can take. 
Condition (|lj) then becomes the familiar 



Pi In dj = (In a) < 



(13) 



Px t+1 {*) = 



Px t+1 {X) = 



P at (a)da / P bt (b)db 

J — OO 

oo 

P Xt (Y)S(X -aY- b)dY , (14) 



^^da r P bt (b)P Xt ( — )db 
a ./__, a 



(15) 



The two pdfs Px t and Px t+1 approach a common sta- 
tionary pdf, P(X), for large t |Q|5|. We are interested 
in the description of the tail of P(X), i.e. for X 3> 6. 
We can then neglect the b term of Px t ((A" — b) /a) in the 
r.h.s. of (jl^). This allows us to simplify ( |l5| ) into 



P(A) 



^Mp^Jda for large A , (16) 
a a 



using P bt (b t )db t — 1. Since ( |l^ ) is linear in P(A), 
the general solution can be written as a sum over a set of 
particular solutions j3l) . These solutions are composed of 
power laws and faster decaying functions (exponentials) . 
The set of power law solutions is obtained by assuming 
the form P(A) - for A > 1. This yields (§) 

determining the exponent fi. 

The inequality (^) and the equation ^ are the corner 
stones of our analysis. We construct and analyze several 
examples whose parameters are constrained by (|5j) and 
we study the solutions of ([?]) and compare them to direct 
numerical simulations. 



This retrieves the regime (||) discussed above. Usually, 
IFS are studied in situations where all the affine transfor- 
mations have their Lipschitz constant individually neg- 
ative, i.e. all are contractive. The present work (where 
D = 1) deals with a rather special but very interest- 
ing situation where some of them are dilating, while on 
average the set of transformations is contractive. This 
correspondence and the discovery that power law distri- 
butions are found when some of the transformations of 
the IFS are dilating suggests to us an investigation of the 
behavior of similar intermittent dilating IFS in higher 
dimensions, where rotations are added to the translation 
and dilation processes. This is left for future work. Here 
we will next use the correspondence with IFS to under- 
stand intuitively the fractal structures found when the at 
take a finite number of values. 



D. Probability density function 

Calling P at , p, f , and Px t+1 the pdfs of at, bt, and 
A t+1 respectively (and assuming that they are integrable 
functions) , then the pdf of X t (obtained by the standard 
Markov argument) obeys the following equation : 



III. TRANSITION OPERATOR APPROACH 

One of the obstacles in the implementation of the ap- 
proach in the previous section is that we need to assume 
that all the considered distributions have densities, which 
is not the case with at least some of our examples. To 
study a more general situation, let us consider the so- 
called transition operator approach. 



A. Transition operator approach and non-smooth 
distributions 

According to our definitions, ([!]) defines a Markov 
chain. We can define the transition operator P of this 
random process as : 



Ph(x) 



h('- 



) dV a ,b{a, b) 



(17) 



for any integrable function h. This operator describes 
the image of a distribution density under the action of 
our random process. If an invariant probability density 
P exists it should satisfy 



4 



pp 



(18) 



The introduction of (|17|) is justified by the fact that this 
integral operator allows for the study of non-smooth and 
even discontinous distributions. 

Consider a simple implementation of a random selec- 
tion scheme. Assume that < <zi < 1 < a%, < b\, 
< 62 for the two maps : 



x — > a\x + bi 
x — > a 2 x + b 2 



(19) 
(20) 



such that ( p9| ) is chosen with the probability p, while 
( ^p| ) is chosen with the probability 1 — p. The corre- 
sponding joint distribution V a ,b for the random variables 
at and b t is singular. Furthermore, these two random 
variables heavily depend on each other. Therefore ( |l5| ) 
cannot directly be used. However it is trivial to specify 
the corresponding transition operator : 

Ph{x) = JL^Szh) + iz£ fc( £z*») . (21) 

di ai a 2 a 2 

Moreover, this description is easily generalized to the case 
with an arbitrary (albeit finite) number of linear maps 
+ bi , where each is selected with the probability 
Pi (provided Y,iPi = !)> 



Ph(x)=J2-H 

— ' n ■ 



Pi. ,x- b,. 



(22) 



The representation fl22|) shows that our random system 
lacks a smooth (and even bounded) invariant density for 
any choice of the positive coefficients ai and bi . Assume 
on the contrary that such pdf h(x) exists. This has to 
satisfy Ph(x) = h(x) for any x G [0, 1]. To proceed fur- 
ther, we need to estimate the variation of the image of 
the h function. The variation of a function is, roughly 
speaking, an integral of the modulus of the derivative of 
the function over its domain. In the case of monotonic 
functions, it can be shown that 



var(P/i) = var 



According to (113) , we note that 



IK <i 



which yields 



> 1 



(23) 



(24) 



(25) 



Therefore, each time we apply the transition operator, 
the variation of the image of a function is multiplied by 
the factor — > 1. Hence the limit distribution (if it 
exists) cannot be a function of bounded variation. 
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FIG. 3. The (numerically obtained) pdf for a t with a two 
point distribution (section at a — 2, ^ = 2, p = 0.95 and 
bt = 1 (10 s iterates, 10 4 equispaced bins per unit of logXt). 

A special case was treated in [^l] with a finite system 
of random maps producing a discrete invariant distribu- 
tion : 



x — * 1 with probability p , 
x — > ax + 1 with probability 



l-p 



(26) 
(27) 



In this case, it is easy to find the invariant distribution an- 
alytically. However, since the first map has a zero value 
of the multiplier a, the system does not satify our as- 
sumptions that all coefficients should be positive. 



B. Existence of pdf for the case of the smooth 
distribution of coefficients 



In this subsection, we study a more general case, where 
we have a random map x — > ax + b with random coeffi- 
cients a, b, whose joint probability distribution is V(a, b). 
The case considered above corresponds to the discrete 
distribution V(a, b). Our main aim here is to prove that, 
if the distribution P(a, b) has a density p(a, b) with "good 
enough" properties, then the random map system also 
has a finite invariant density. Indeed, consider the corre- 
sponding transition operator: 



Ph(x) 



p{a,b) x 



- ) dadb 



(28) 



After the change of variables £ = x — b this operator 
becomes 



Ph(x) 



Assume now that 



Pia > x -® h(i)dad£. (29) 
a a 



var(,p(a, b) < C < 00 



(30) 
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for any a > 0, where v&Tbp(a, b) stands the variation of 
the p(a, b) function with respect to the second variable. 
Using the above representation, we find that 



var(P/i) < J J ^-h{^)dadi = C J 



C( / h(x) dx) I da . 



h(-)d(~) | d„ 
a a 

(31) 



Since h is assumed to be the density of a probability dis- 
tribution, then J h(x) dx is finite. This shows that the 
variation of P/i is universally bounded from above. The 
existence of the invariant distribution |33| ] then proves 
the existence of the pdf. 



C. Markov dependent choice of the subsequent map 

Our earlier discussion of random map systems assumed 
that the choice of a subsequent map x — > aiX+bi does not 
depend on the immediatly antecedent map chosen. This 
is a relatively strong restriction and in this section we 
shall show that this assumption is not necessary. We will 
show that, for the stationary process, the representation 
of the transition operator depends only on the station- 
ary probabilities of the random choices pi and not on the 
transition probabilities between subsequent maps. 

Assume that currently the map x — > a L x + bi was cho- 
sen, then the conditional probability to choose the map 
x — > cijx + bj is equal to Pij. The process of the ran- 
dom choice is governed by the finite state Markov chain 



with the transition probabilities (pij). Assume that this 
Markov chain is ergodic and denote by pi its unique in- 
variant distribution. Then our entire system is still a 
Markov chain, whose transition operator is: 

= j2-H x —^)J2p.P,, = '£- h i—'i^ < 32 > 

^ — ' n j n ■ — ' * — ' n j a _• 



It depends only on the stationary probabilities pi . As a 
result, we immediately see that all asymptotic properties 
also depend only on the choice of pi . 

The generalization of our argument for the general case 
where the joint distribution of the coefficients a and b 
may have both discrete and continuous components is 
straightforward. 



IV. TWO-POINT DISTRIBUTIONS 

Let us return now to the question about the asymp- 
totic (as x — > oo) properties of pdf for the case of only 
two linear maps. The above derivations show that these 
asymptotic properties do not depend on the choice of 



the additional terms bi (as long as they are non-zero and 
well-behaved). Let 



a > 1, <p < 1, and £ > 



(33) 



such that a\ = 1/a < 1 and a-z = cfi > 1. Equation (|l6| ) 
becomes 

P(X) = paP{aX) + (1 - p)a,-tp(a- ( X) . (34) 
Condition (§) imposes the requirement 



and equation (0) leads to 
(l-p)z 1+ « - z+p = 



<P< 1 , 



(35) 



where z = o/ 1 and z complex . (36) 



A. £ integer 



At first take £ = 1, then from fl3q), we see that p must 
be within i < p < 1. The two real solutions of (|36j) are 

z± = where A = 1 - 4p(l - p) > 0. From the 

definition z = a M (and e l2n7r = 1 for any integer n), we 
obtain 



In z± 2im 

M±,n = M-R + W = -; 1" i--. 

in a ma 

The pdf of X t is thus of the form 



a 



±,n 



±,n 



X 



cos(^/ In X t 



(37) 



(38) 



The prefered scaling ratios are obtained by the factors of 
X t reproducing the same values of the cosine, i.e. a~ , 
with n integer. The discrete scale invariance is simply 
the result of the intermittent amplification by the fixed 
factor a. The log-periodicity is thus trivially associated 
with the discrete multiplicative structure. 

When £ = 2, three real z solutions exist for z, for the 
allowed range | < p < 1. The imaginary part of [i thus 
stems from the same technical reason as for £ = 1 and 
reflects the intermittent amplification by the factor a 2 . 

In general if £ = N is a positive integer, (^) obeys 

(l-p)z N+1 - z + p = Q & {z-l)Q N (z;p) = , 

(39) 

where 



N 
\ " k 



Qn(z;p) = (l-p) 2_^z 
fe=i 



P ■ 



(40) 



G 



The root structure of Qn(z;p) — for SKz > is deter- 
mined with e.g. Routh's algorithm (4^]. It can be shown 
that for N < 5, this polynomial has only one root with 
dlz > and that 3z = for this single root. However 
for N > 5 there always exist roots such that 5Rz > and 
3z ^ 0. 

To illustrate the integer £ regime, Fig. [l] shows a seg- 
ment of the AT t history for the case a = 2, £ = 2, p = 0.95, 
and a constant bt = 1. Most iterates are small while rare 
intermittent excursions explore very large values. The 
(numerically obtained) cumulative distribution is shown 
in the log-log plot of Fig. 0. 
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FIG. 4. The scaled pdf of X given by (gj) for a t with two 
point distribution at a — 2, £ = 2, p = 0.95, and b t uniform 
with (a) 1 x pdf and (3 = || (10 8 iterates, 10 4 equispaced bins 
per unit of \ogX t ), (b) 10^ x pdf and /3 = | (10 9 iterates, 10 3 
equispaced bins per unit of logXt, the same for (c-d)), (c) 
10 6 x pdf and = §, (d) 10 9 x pdf and /3 = 0. 



A complex structure, reminiscent of a devil staircase, 
overlays an average linear decay. The structure corre- 
sponds to all possible values of n in the imaginary part 
( |37j ) of the exponent [i, where the largest n provide the 
smallest details of the cumulative distribution. Figure || 
shows the (numerically obtained) pdf, i.e. the deriva- 
tive of Fig. g. We observe a self-similar structure, as ex- 
pected fro m th e correspondence with the IFS discussed 
in section |II C| (IFS in general encode stochastic fractal 
structures J45[| ) . We also observe that the distribution 
seems to be nowhere continuous, as expected from the 



It is interesting to progressively coarse-grain this self- 
similar structure by introducing a disorder on 6 t . This 
is accomplished by chosing b t uniformly in the inter- 
val \J3, 1] . The value j3 = 1 recovers the ordered case 
b t = 1. Decreasing (3 corresponds to increasing the dis- 
order. Figure [I| shows the pdf P(Xt) for decreasing 
values (3 = ||, |, | ,0 (while keeping a — 2, £ = 2, and 

p = 0.95). The roots of (f|) are z = 1 and z± = ~ 1± 2 v/T7 , 
or v R0 = 0, ,U/o = n-jj^, n R+ » 1.9588, n I+ = Hio, 
fiB~ ~ 2.5890, jti/_ = (1 + 2n)j^2 for integer n. 




frequency 

FIG. 5. (a) The logarithmic derivative of a portion of the 
pdf tail for the upper (/3 = 0) trace in (^). (b) The Lomb 
periodogram of (a). 



by 



The "frequency" of a log-periodic oscillation is defined 



/(") 



1 





27r log e In A 



(41) 



where we define A as the scaling ratio associated with 



the log-periodicity |LJ2J. Here /(1) = /(1)+ « 3.3219, 
/(0)_ w 1.6609, and /(1)_ w 4.9828. These numbers are 
compared to the spectrum analysis of the tail portion of 
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the pdf. We use the logarithmic derivative of the [3 — 
pdf of Fig. |H to get a data with zero average slope (its 
averave value is the leading power law exponent). This is 
represented in Fig. |^ a. Its Lomb periodogram spectrum 
is shown in Fig. || b and yields four frequencies 0.6, 
1.6, 3.3, and 5.0. The smallest of these is attributed to 
the inverse of the (log) tail length used. The others are 
in good agreement with the predictions /(0)_, /(1)+, 
/(l)- in ascending order. The small, partially hidden, 
bump at 2.1 and the more recognizable one at 4.3 are of 
unknown origin. 



B. £ non-integer 



10 L 
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10000 



FIG. 6. (a) The pdf for at with two point distribution at 
a = 2, £ = 2.5, p = 0.95, and b t = 1 (10 8 iterates, 10 4 equis- 
paced bins per unit of logXt). (b) The scaled pdfs for at with 
two point distribution at a — 2, £ = 2.5, p = 0.95, /3 = jg 



(lower trace, 1 x pdf), and /3 = (upper trace, 10 3 x pdf). 
Both > 
logXt. 



Both with 10 9 iterates and 10 3 equispaced bins per unit of 



We have shown that for integer £ > 5 there will al- 
ways exist some roots of (|3rj ) with non-zero imaginary 
and positive real parts. This type of root structure is 
common for non-integer £. If £ is irrational then an infi- 
nite number of distinct roots solves ([56j). We select the 
slightly simpler case £ = 2.5 (with a = 2 and p = 0.95 as 
before) . 

In Fig. H a, the pdf for bt = 1 is given for the 10 8 
first iterates of ([!]) with a binning density of 10 4 points 
per decade. Figure ^ b shows the pdf for a uniformly 
distributed b t with (3 = y| (lower trace) and the scaled 
pdf (upper trace, 10 3 x pdf) for (3 = 0, i.e. b t uniformly 
distributed between and 1. These two pdfs use the first 
10 9 iterates of (|l|) with a log-equidistant binning of 10 3 

points per decade. 

As already pointed out in section [II B , a continuous 
bt distribution seems to lead to a continuous pdf for X. 
The analysis above neglected the influence of a varying 
b t . We resort to a limit consideration on a sequence of 
progressively thinned, uniform, bt distributions to match 
theory with the present simulation results. We select the 
three cases (3 = 0, |, |. For each of these tail regions of 
the pdf (Fig. a, Fig. H a, and Fig. ^| a), the logarithmic 
derivative is computed (Fig. [?| b, Fig. || b, and Fig. ^| b). 
This gives a local estimate of the leading exponent of 
the power law tail of the pdf. A constant value would 
correspond to a pure power law. Oscillations which are 
approximately periodic in In AT are the signatures of the 
log-periodicity. This is confirmed by a_spectral analysis 
given in Fig. Mc, Fig. 
shown in Fig.ll] b, Fig. 



c, and Fig. He of the signals 
I b, and Fig!f|_b respectively, 
using the Lomb periodogram technique |46|| . We clearly 
identify a number of frequencies. We compare these nu- 
merical results with a direct analytical determination of 
the roots of (|3|3^). The complex \i solutions are sought 
where z = a M = e x+iy . These are the roots of 



(42) 



where 



R(x,y;t,p) = (1 -p)e^ x co S ((l + Qy) 

- e x cos y + p , (43) 

J(x,y;Z,p) = (1 -p)e^ x sm((l +$y) 

-e x smy. (44) 



The zeroes of ([43|) and (J44J) define nodal curves in the x-y 
plane. The solutions (x, y) are the intersections of these 
nodal curves. When £ is rational £ = M/N, where M and 
N are the smallest relative prime positive integers, we see 
that the set of solutions is periodic in the ^/-direction with 
a period 27rA. The two equations (^,Q) are converted 
into 



(l-p)e^ 



sin y 



e = p- 



sin(l + £)y ' 
sin(l + Qy 



(45) 
(46) 
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This shows that the values of x are bounded from above 
by a finite number. An unbounded x would from ( Rj ) 
correspond to a vanishing sin£y. This would imply 
that £y — nir, for some integer n, and therefore that 
sin(l + £)y = (— l) n sint/. Using ( ff5| ) we would obtain 
(1 — p)e^ x = 1 (n odd is not allowed) leading to a con- 
tradiction. This implies that there is a maximum value 
for \xr. Table ffl gives the five solutions (x m , ym)i indexed 
by m = 1 to 5. For a given m, we also extract the few 
first solutions obtained by An periodic repetitions in the 
y direction. These solutions are indexed by an additional 
integer n corresponding to the order of the 47r period. We 
give these solutions in the /ir and /i/ parameter space. 

TABLE I. All roots x, y, and the first few [i = (ir + 
roots (here a M = e x + l vj Q f Q for a t with two point 
distribution at a = 2, £ = 2.5, and p = 0.95. 



m 




1 


2 


3 


4 


5 


x(m) 




1.03 


1.28 


1.19 


1.19 


1.28 


y(m) 




0.0 


2.56 


4.91 


7.66 


10.06 






1.47 


1.85 


1.72 


1.72 


1.85 


fj,i(m,n = 


0) 


0.0 


3.69 


7.08 


11.04 


14.43 


fj,i(m, n — 


1) 


18.13 


21.82 


25.21 


29.18 


32.57 


[ii(m, n — 


2) 


36.26 


39.95 


45.34 


47.31 


50.69 


fii(m, n — 


3) 


54.39 


58.08 


61.47 


65.43 


68.82 



The appendix provides an approximate but quite ac- 
curate analytical determination of the solutions found in 
this table, based on a perturbative scheme. The pdf of 
X t is thus a sum of power laws overlayed by log-periodic 
oscillations of the type 



P(Xt) = 



X 



cos(/x/(m, n) \n.X t ) . 



(47) 



The leading power law behavior is given by the first m = 
1 real solution, which has the smallest [i = /ii?.(l) ~ 1-47 
(with £t/(l,0) = 0). The other solutions have larger 
and thus corresponds to subleading corrections. We de- 
fine the "gap" as the smallest difference between the real 
parts of the complex solutions to the first real solution. 
This gap measures the strength of the subleading log- 
periodic corrections to the leading power law behavior. 
In the present situation, all fiR(m) take two values ~ 1.72 
and « 1.85, which are close to /ir(1). The gap is approx- 
imately 0.25. This corresponds to strong corrections to 
the leading scaling for which the log-periodic oscillations 
are very visible. Notice that, asymptotically, the oscil- 
lations disappear for X t — > oo, as /x«(m > 1) > /Ur(1). 
This effect is very weak in the present case, since the rel- 
ative amplitude of the dominant log-periodic oscillations 
decays as X^ 



"D 
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logX, 




4 

frequency 

FIG. 7. The pdf tail for at with two point distribution at 
a = 2, £ = 2.5, p = 0.95, and b t uniform with (3 = (10 9 
iterates, 10 3 equispaced bins per unit of logXi) (a), (b) Its 
logarithmic derivative, (c) The Lomb periodogram of (b). 
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The pdf tail for at with two point distribution at 

n9 



FIG. 

a = 2, £ = 2.5, p = 0.95, and b t uniform with /3 = § (10 M 
iterates, 10 3 equispaced bins per unit of \ogX t ) (a), (b) Its 
logarithmic derivative, (c) The Lomb periodogram of (b). 
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FIG. 9. The pdf tail for at with two point distribution at 



2, £ = 2.5, p = 0.95, and b t uniform with f3 



(10 y 



iterates, 10 3 equispaced bins per unit of logXt) (a), (b) Its 
logarithmic derivative, (c) The Lomb periodogram of (b). 
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TABLE II. The predicted and the observed frequencies 
(obtained by spectral analysis of Fig. |?] c, Fig. ^ c, and 
Fig. c). The increased disorder on the variable bt is noted 
with the different (3 subscript. The first and the nth rows 
(the F(M) labelled rows) contain the predicted frequencies. 
The F-f3 subscripted rows list the frequencies retrieved from 
the numerical realizations. The bold type is used to indicate 
the frequency that gives the largest peak in the spectrum 
window, whereas the other well-defined peaks are given in 
standard type. The parenthesis enclosed values correspond 
to spectrum peaks that barely are above the noise level of the 
Lomb periodgram. 



where 



M 


1 


2 


3 


4 


5 


F(M) 


1.35 


2.60 


4.05 


5.29 


6.64 


Fg= r /a(M) 






4.13 


5.25 


6.63 


F 0=3 /4 (AO 


(1.33) 


(2.53) 


4.08 


5.25 


6.65 


Fp =0 (M) 


1.31 


2.64 


4.08 




(6.69) 


M 


6 


7 


8 


9 


10 


F(M) 


8.00 


9.24 


10.69 


11.93 


13.29 


F 0= r/ S (M) 


(7.83) 


9.13 




11.96 


13.25 


*>=3/4(Af) 


(7.84) 


9.12 


10.77 


11.92 


13.30 


ff»=o(M) 













Table || gives the observed frequencies obtained by the 
spectral analysis of Fig. c, Fig. [sj c, and Fig. ^| c and 
compare them with the predicted values. This contains 
the different cases with increased disorder on the vari- 
able bt- Here a slight generalization of frequency is used, 
f(m,n) ee and F(M) = f(M+l-5n,n= [f ]). 

Increasing the disorder in the variable bt increases the 
noise level and progressively washes out the higher fre- 
quencies. 



V. TWO LEVEL "STAIRCASE" DISTRIBUTION 

We now study a situation with a much larger disorder 
where the multiplicative factors at are selected from a 
broad, continuous, distribution. To minimize the num- 
ber of control parameters for the pdf's, we use distribu- 
tions which are constant by parts. In the next section we 
will examine the uniform distribution. Here, we divide 
the interval a^] into two sub-intervals 1] and [1, a^] 
having different weights p and 1 —p respectively. The idea 
is to allow for a different weight of the damping versus 
amplificating processes and examine the consequence on 
the amplitude of the log-periodic structures. This choice 
corresponds to the following pdf for at : 



P at (at) = 



P 



1-1/ a 
a«-l 



where O is the Heaviside function 
dition (||) that (In at) < reads 

P(a,Q <P < 1 , 



(6(ot - l/o) - G(at - 1)) 

(0(ot - 1) - Q(a t - a«)) , (48) 
The stationarity con- 



(49) 



with 



z«(£lna- 1) + 1 



a«(£Ina- 1) + l + 5(a,£) 
a« - 1 

8(a,0 = (a-l-lna). (50) 



a - 1 

The integral equation (|lj]) is now 



P(X) = 



ap 



a-lJx/a a t 



1 



da* ■ 



(51) 



This equation has a power law solution for large X if the 
exponent fj, satisfies (0), leading to 



M+ 1 



ap 
a-l 



(1 - a 



1-P 
a«-l 



(d 



1) 



(52) 



Assuming a complex solution /i = (ir + ifij splits equa- 
tion (B2j) into 



ap 
a- 1 
1-p 



+ 



ap 



I _ a -(MR+l) C0S ( M/ mfl )) 

( a £0**+i) cos(£/// In a) - 



1 

-(mr+i) 



a - 1 
1-P 



sin(/i/ In a) 



-a 



1 



sin(£/x/ In a) 



(53) 



(54) 



To allow for a comparison with the previous case, we 
keep the same parameters a = 2, £ = 2.5 and p — 0.95 
as before. The solutions of these equations are graphi- 
cally represented as the intersections of the continuous 
and dashed lines in Fig. [Ic] (here p(2.0,2.5) w 0.78266). 
Table III lists the smallest roots and their corresponding 
log-periodic frequencies (f(m) = Ml 2 ( -"- > ) . 

20 



=E 10 




FIG. 10. A portion of the complex /x plane with roots 
p, = fiR + ifi[ of (|5^) (for a t with a two level staircase distri- 
bution at a = 2, £ — 2.5, p — 0.95) as intersections between 
its real part (solid lines) and imaginary part (dashed lines). 
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10000 



FIG. 11. The pdf for the case where the pdf of at has a 
level staircase structure at a = 2, £ = 2.5, p = 0.95, with 
b t — 1 (10 9 iterates, 10 3 equispaced bins per unit of logXt). 



TABLE III. The first few fi = fj, R + ifc roots of (|5J, |54|) 
and the predicted log-periodic frequencies for a t with a two 
level staircase distribution at a — 2, £ = 2.5, and p = 0.95. 



m 


1 


2 


3 


4 


5 


6 


p,i(m) 
f(m) 


1.6535 
0.0000 
0.0000 


2.4918 
4.3475 
1.5932 


2.8418 
8.0152 
2.9373 


3.0389 
11.6794 
4.2801 


3.2015 
15.3255 
5.6163 


3.3179 
18.9598 
6.9482 



An important difference with the previous two-point 
pdf is that now the gap value is ^r(2) — (j.r(1) — 0.84, 
which is about three times larger than before. This 
means that the log-periodic structure are smaller and 
decay faster for large X. They are still quite visible as 
found in Fig. [Il], where we can observe the undulation 
of P(X) , s tail. The results obtained for the various bt 
distributions (from a non-random bt — 1 to a uniform bt 
with (3 = 0) are essentially the same for X > 3. The 
only difference is that a larger disorder in b t allows for an 
exploration of the interval closer to 0. The Lomb power 
spectrum analysis is presented in Fig. [l|. The fundamen- 
tal frequency, /(2) w 1.59, is visible in all the simulations 
and more clearly in the spectral analysis where a strong 
peak appears in the Lomb power spectrum. The next 
higher frequency, /(3) ~ 3.0, is the only one which can 
be detected as the disorder in b t increases. All higher 
frequencies are lost in the noise. The reason for this is 
clear. The relative amplitude of a given frequency f(m) 
is quantified by /i^(m) — /Lt^(l). For the second frequency, 
we have [Ir(3) — /4r(1) = 1-19. For the third frequency, 
we have ;U.r(4) — /J,r(1) = 1.39 and so on. It seems that 
a difference fiR{m) — (J,r(1) of the order or less than 1 
is necessary for the clear detection of log-periodicity. In- 
tuitively, this ensures that the amplitude of the oscilla- 
tions does not decay more than by a factor 100 over two 
decades. We notice that a similar gap about 1 was found 
in the analysis of the log-periodic structure of DLA clus- 
ters |17) . The present analysis rationalizes why we have 



been able to detect these structures in this case. 
-2.5 
-3.0 




2 3 
frequency 

FIG. 12. (a) The Pdf tail for at with two point distribu- 
tion at a = 2, ^ = 2.5, p = 0.95, and b t = 1 (10 9 iterates, 
10 3 equispaced bins per unit of logXt). (b) The logarithmic 
derivative of (a), (c) The Lomb periodogram of (b). 



Table IV makes the comparison between the predicted 



and observed frequencies. 
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TABLE IV. Predicted and observed frequencies for a t with 
a two level staircase distribution at a = 2, £ = 2.5, p = 0.95, 
and two different choices for bt distribution. The first row con- 
tains the predicted frequencies. Subsequent rows are the fre- 
quencies retrieved from the numerical simulations. The bold 
emphasis indicates the frequency that gives the largest peak 
in the spectrum. The other well-defined peaks are written in 
normal format, while numbers inside parenthesis correspond 
to peaks in the spectrum that are barely above the noise level 
of the Lomb periodgram. 



m 


1 


2 


3 


/(m + 1) 


1.59 


2.94 


4.28 




1.51 


3.07 


(4.31) 


fp=o(m) 


1.51 


3.00 


(4.50) 



We conclude that a log-periodic structure of the tail 
of the X t 's pdf is present for the smeared out two level 
staircase distribution, although its amplitude is weakened 
as compared to the previous two point distribution case. 
This was expected from the theoretical analysis of the in- 
fluence of disorder (2^jlJ]. The important aspect of our 
result is that the log-periodicity, and the prefered scal- 
ing ratios A, can no more be associated to a specifically 
chosen amplification factor, as for the previous two-point 
distribution. Notwithstanding the presence of a large 
disorder, a discrete set of effective scaling factors are se- 
lected. It is amazing to us how strong is this effect and 
how relatively weak is the influence of the disorder. 



VI. UNIFORM DISTRIBUTION 



This last section deals with the effects of a very strong 
disorder on a t . To attain this goal, we consider a uniform 
a t -distribution 



9(a t - at) - 6(a t - a r )) 



a r — at 



(55) 



where 



< a e < 1 , 

1 < a r < a r (ae) 



(56) 
(57) 



a r (ae) is such that (||) is obeyed and is the solution of 

a r (ae) In a r (ae) — a r (ae) = a^lnci£ — ae . (58) 




FIG. 13. Allowed ai~a r domain for at with uniform distri- 
bution. The dotted boundary depicts the strict inequalities 
given in (|5t|) . 

The largest width compatible with (^) corresponds to 
ae = and a r = e. When ae is made larger, the maxi- 
mum value of a r progressively decays towards 1 . 

The integral equation (|l^) is 



P(X) 



a T p(2L) 
, a* 



(59) 



The tail of P(X) takes the form of a power law if the 
exponent [i is the solution of 



(ji - l)(a r - a t 
With /i = /Xfl + we get 

((1r + l)(a r — ae) = a% R+ cos(/j,/ lna r ) 

— a^ R+1 cos(/z/ In ae) 
Hi(a r — ae) = a^ R+1 sm(/ij lna r ) 

— a^ R+1 sin(/i/ In a^) . 



(60) 

(61) 
(62) 



We select ae = 0.001 and a r = 1.9. The solutions with 
the smallest, positive, real parts are given in the table [v|, 
together with their corresponding log-frequencies /. 

TABLE V. The first few fj, = fm + ip,x roots of @ ||) 
and the predicted log-periodic frequencies for the uniformly 
distributed a t with ae — 0.001 and a r — 1.9. 



m 


1 


2 


3 


4 


5 


fii(m) 
f(m) 


1.2667 
0.0000 
0.0000 


3.9414 
11.6095 
4.2545 


4.8349 
21.6147 
7.9211 


5.3984 
31.5024 
11.5446 


5.8116 
41.3494 
15.1532 



Figure 13 shows the allowed ae~a r region. 



The most striking feature to note is the large gap value, 
(J,r(2) — jUfl(l) — 2.67. These differences increase rapidly 
with the order m of the solution. This implies that the 
oscillations must be extremely weak and severely damp- 
ened. In Fig. |lj we explore the dependence of the gap as 
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a function of the parameters at and a r of the model. If 
a£ — and a r = 2.71, we find /Ur(2) — /ifl(l) — 2.1, but 
/ifi(l) = 0.006 is very small. When = and a r = 2.0, 
we find /i fl (2) - /i fl (l) = 2.545 with = 1.0. 




FIG. 14. The gap value, ^r(2) — fj,n(l), as a function of 
a r for different choices of a;. 



The situation does not improve if we take at — > 1~ and 
a r — > 1+ (while keeping the stationarity condition (||)). 
Notice that we need a r > 1 in order to get a solution for 
/i, i.e to get a power law pdf for X t . This stems from the 
fundamental fact that the power law pdf results from in- 
termittent amplifications. In summary, the log-periodic 
oscillations are present theoretically but are very difficult 
to measure and quantify. The pdf for the bt = 1 case is 
shown in Fig. |l5| . 




10000 



FIG. 15. The pdf for at with uniform distribution between 
ai = 0.001 and a r = 1.9; bt — 1 (10 9 iterates, 10 3 equispaced 
bins per unit of log Xt). 

A segment of its tail is analyzed with the same proce- 
dure as was applied to the previous two at-distribution 
families (Fig. 16). 




10 

frequency 

FIG. 16. (a) The logarithmic derivative of the pdf of X 
for at with a uniform distribution between a; = 0.001 and 
a r = 1. 
logXt) 



, b t = 1 (10 iterates, 10 equispaced bins per unit of 
(b) The Lomb periodogram of (a). 



Only /(4) seems to be recognizable among all the peaks 
in the Lomb periodogram but this seems even far-fetched 
as the signal is within the noise level. In conclusion, 
this uniform case corresponds to a large gap and the log- 
periodic structures, that are present in theory, are not 
clearly visible. 



VII. CONCLUDING REMARKS 

We started our analysis by considering the intermittent 
multiplicative processes of a simple binomial at distri- 
bution. Not surprisingly, strong log-periodic corrections 
to the main power law probability density function have 
been found for the random affine map. The disorder in 
the additive constant smooths out the higher frequencies 
but does not dampen out the smallest log-periodic fre- 
quencies. We then analyzed situations with increasing 
disorder in the multiplicative terms. Instead of going to 
the weak disorder regime with two broadened peaks, we 
analyzed a pdf of multiplicative factors which consists 
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of a two-step staircase. In this already large disorder 
regime, we have found that the log-periodic structure of 
the tail of the X t 's pdf is present for the smeared out two 
level staircase distribution, although it is weakened in its 
amplitude as compared to the two point distribution case. 
The most important aspect of our results is that the log- 
periodicity, and the prefered scaling ratios A, can no more 
be associated to a specifically chosen amplification fac- 
tor, as for the two-point distribution. Notwithstanding 
the presence of a large disorder, a discrete set of effec- 
tive scaling factors are selected. The "gap", defined as 
the difference between the smallest exponent real part 
and the real solution, controls the strength of the log- 
periodicity. We have been able to determine that the gap 
must be of the order or less than 1 in order for the log- 
periodicity to be strong. Larger gaps still lead to visible 
effects but the analysis must then be very precise and the 
noise level very low. This is the situation found for a uni- 
form distribution of multiplicative factors. In summary, 
we have shown that log-periodicity remains a significant 
effect even in the presence of significant disorder. 



VIII. APPENDIX : DETERMINATION OF THE 
EXPONENTS FOR THE TWO-POINT 
DISTRIBUTION USING A PERTURB ATIVE 
ANALYSIS 

The power law structure of the P(X) pdf character- 
izes rare excursions of X t to large values. These large 
values are reached by repeated occurrence of the ampli- 
fying multiplicative factor . This motivates us to make 
the approximation of neglecting the first "damping" term 
in the r.h.s. of (|34|), 



P(X) w (1 -p)a~^P{a~ i X) . 



(63) 



This functional equation is simpler to handle. It also has 
a form reminiscent of the renormalization group equation 
that Feigenbaum used in his analysis of a bifurcation se- 
quence of the logistic equation W7| . Analoguous equa- 
tions has also been discussed in |1 l| 4jg|]2l| . Assuming a 
power law form for P{X) (P(X) = x u^ ) provides for 
the y, equation : 



(1 -p)a^ = 1 . 
With the notation z — e x+vy — a MH+w , we obtain 



(64) 



./■ = — — ln(l — p) and (j,r = — — — (65) 



2kn 2?r 
V = —r~ and \xi = 



£ lna 



t hiAfe 
where A& = cfi' k , k integer 



(66) 



An approximation of the imaginary part of the roots of 
([l2|) therefore is 



2fc7T 



(67) 



This agrees fairly well with the exact computed roots 
shown in table || We can improve on this estimation by 
inserting the parametrization 



y 



2klT 



(68) 



in the full set of equations ( J4q ) and (f44|). This yields the 
following system of two equations of the two unknown x 
and e : 



((1 - p)eS 1+ ^ x cos e(l +0~e x cos e) cos 
-((l-p) e ( 1+ « )a; sme(l + £) 



2kir 



2/C7T 

e x sin e) sin — - — h p = , 



(69) 



((1 - p)e^ x cos e(l + £) — cos e) sin 



2fc?r 



+ ((l-p)e t *sine(l + £) 
2kn 

— sine) cos — — = . 



(70) 



Assuming e(l +£) to be "small", we expand the trigono- 
metric functions to first order in e : 

((1 - p)e {1+ t )x - e x ) cos2fc7r/e - e((l - p)(l + i)e {1+ ^ )x 
-e a; )sin2fc7r/^ + p = 0, (71) 

((1 ~p)e^ x - l)sin2A;7r/C + e((l-p)(l + 0e Sa! 

- l)cos2fc7r/£ = . (72) 



Eliminating e between the two preceding equations, we 
get an equation in the sole variable x : 



(1 -p)e {1+ ^ x - e x +pcos2/br/f; = 



(73) 



There is a unique solution in x for each k. Knowing x, 
we then get e from 



(l-p)e« 



1 



2kn 

(i-p)(i + 0^-i tan X- (74) 



Table VI compares these solutions to the exact ones, in 
the case when a — 2, £ = 2.5, p = 0.95. 

TABLE VI. The appoximate (according to (fr|, |^, |§)) 
and the exact (from fil|)) roots for a t with two point 
distribution at a = 2, £ = 2.5, and p = 0.95. 



k 







x ap pr X {k) 1.0333 1.2760 1.1597 1.1597 1.2760 1.0333 

Xexact(k) 1.0333 1.2808 1.1922 1.1922 1.2808 1.0333 

e(fc) 0.0000 0.0480 -0.1301 0.1301 -0.0480 0.0000 

y apP rx(k) 0.0000 2.5612 4.8965 7.6699 10.0051 12.5664 

yexact{k) 0.0000 2.5605 4.9101 7.6563 10.0058 12.5664 
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